Vortex arrays in type-II superconductors admit the translational symmetry of an infinite system. There are cases, however, like ultra-cold trapped Fermi gases and the crust of neutron stars, where finite-size effects make it quite more complex to account for the geometrical arrangement of vortices. Here, we self-consistently generate these arrays of vortices at zero and finite temperature through a microscopic description of the non-homogeneous superfluid based on a differential equation for the local order parameter, obtained by coarse graining the Bogoliubov-de Gennes (BdG) equations. In this way, the strength of the inter-particle interaction is varied along the BCS-BEC crossover, from largely overlapping Cooper pairs in the BCS limit to dilute composite bosons in the BEC limit. Detailed comparison with two landmark experiments on ultra-cold Fermi gases, aimed at revealing the presence of the superfluid phase, brings out several features that makes them relevant for other systems in nature as well.